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Section  1:  Need 

The  new  generation  of  fighter  aircraft  include  internally  carried  stores.  Fighter  aircraft  such 
as  the  F-22  and  the  Joint  Strike  Fighter  must  separate  weapons  over  a  wide  range  of  Mach  num¬ 
bers  and  maneuver  conditions.  The  acoustic  and  aerodynamic  loads  for  low  altitude,  supersonic 
weapons  release  are  close  to  an  order  of  magnitude  higher  than  is  found  for  high  altitude,  tran¬ 
sonic  weapons  release.  As  a  result,  predictions  of  acoustic  loads  and  store  separation  weapons 
bavs  has  taken  on  new  importance.  The  cost  of  wind  tunnel  testing  for  assessment  of  acoustic 
loads  and  store  separation  in  the  development  phase  of  a  fighter  aircraft  program  can  run  into  tens 
of  millions  of  dollars.  Subsequent  flight  testing  for  certification  requires  additional  tens  of  mil¬ 
lions  of  dollars.  Development  of  computational  methods  for  acoustic  analysis  and  weapons  sepa¬ 
ration  prediction  holds  the  promise  of  saving  many  millions  of  dollars  in  wind  tunnel  and  flight 
testing  and  reducing  design  cycle  time  on  future  aircraft  programs. 

The  internal  weapons  bay  of  a  modem  fighter  aircraft  is  extremely  complex  geometrically 
due  to  the  pre.sence  of  bulkheads,  stores,  bay  doors  and  the  storage  of  a  variety  of  flight  hardware 
in  the  bay.  The  flow  in  and  near  the  weapons  bay  is  highly  unsteady.  The  flowfield  external  to  the 
bay  is  strongly  influenced  by  both  the  internal  and  external  aircraft  geometry.  As  a  result,  simula¬ 
tion  of  the  weapons  bay  flowfield  using  Computational  Fluid  Dynamics  (CFD)  requires  large 
computational  meshes  in  order  to  adequately  resolve  highly  complex  fighter  aircraft  geometries 
and  flowfields.  The  goal  for  this  AFOSR  research  program  is  to  develop  computational  methods 
that  van  be  used  efficiently  in  the  accurate  prediction  of  acoustic  loads  in  weapons  bays  and  store 
release. 


1.1  Objectives 

The  program  objectives  were  to  develop  and  demonstrate  methods  for  the  accurate  simula¬ 
tion  of  weapons  bay  acoustic  loads  and  to  develop  and  demonstrate  methods  for  weapons  separa¬ 
tion  from  a  fighter  aircraft  weapons  bay.  The  acoustic  load  test  case  is  a  Mach  1.2  Weapons  in 
Cavity  (WICS)  study.  The  store  release  simulation  demonstration  is  for  an  AIM-9  missile  release 
from  an  F-22  at  Mach  1.6.  The  methods  developed  for  this  program  and  the  demonstration  cases 
form  a  foundation  for  the  development  of  engineering  design  tools  for  the  accurate  simulation  of 
the  weapons  bay  environment  and  store  release. 
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These  objectives  proved  challenging  given  the  program  budget.  All  of  the  key  goals  for 
methodology  development  were  met.  Acoustic  load  and  store  separation  demonstration  cases 
were  run,  and  preliminary  results  were  obtained.  The  demonstration  cases  identified  parts  of  the 
methodology  that  need  to  be  streamlined  in  order  to  achieve  more  rapid  turnaround  of  the  solu¬ 
tions. 


For  acoustic  predict!^  iS,  the  LES  methods  implemented  in  the  parallel  version  of  Splitflow 
were  tested  on  a  simple  cavity  flowfield.  Numerical  diffusion  effects  and  numerical  flux  schemes 
were  developed  and  investigated.  A  method  was  developed  for  simulation  of  an  essentially  equi¬ 
librium  compressible  turbulent  boundary  layer  for  input  to  an  LES  cavity  solution. 

To  enable  store  separation  simulation,  important  capabilities  were  added  to  Splitflow.  Mov¬ 
ing  grid  terms  were  added  for  grid  translation  and  rotation.  Methods  for  initialization  of  computa¬ 
tional  cells  originating  at  store  boundaries  due  to  the  motion  of  the  store  were  developed.  Six 
degree  of  freedom  dynamics  (6-DOF)  simulation  capability  was  added  to  Splitflow.  The  output 
from  the  6-DOF  routines  was  used  to  specify  the  motion  of  the  store  through  the  flowfield  at  each 
iteration.  An  interface  method  was  developed  to  allow  regions  of  the  flowfield  with  a  body  fitted 
grid  to  interface  with  a  cut  grid. 

1.2  Overview  of  Weapons  Bay  and  Store  Release  Flows 

The  flow  in  the  shear  layer  over  and  in  an  open  weapons  bay  is  highly  unsteady.  Much  of 
the  large  scale  unsteadiness  in  the  weapons  bay  is  triggered  by  acoustic  feedback.  Deflections  in 
the  separated  shear  layer  over  the  top  of  the  bay  generate  pressure  disturbances.  These  distur¬ 
bances  generate  acoustic  waves  that  travel  upstream  through  the  slow  moving  flow  in  the  weapons 
bay.  When  the  upstream  running  pressure  waves  reach  the  front  wall  of  the  bay,  they  reflect  off  of 
the  front  wall  of  the  bay  and  disturb  the  separated  shear  layer  near  the  separation  point.  The  shear 
layer  disturbance  grows  as  it  is  convected  downstream,  resulting  in  the  generation  of  new  pressure 
disturbances.  This  mechanism  can  result  in  a  resonant  feedback  in  weapons  bays  and  in  simple 
cavities.  In  addition  to  the  large  scale  unsteadiness,  the  incoming  boundary  layer  and  separated 
shear  layer  are  highly  turbulent,  include  small  scale  turbulent  structures. 
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Successful  simulation  of  the  weapons  bay  flowfield  requires  capturing  the  effects  of  both 
the  large  scale  unsteadiness  and  the  turbulent  boundary  and  shear  layers.'  Attempts  to  use 
unsteady  Reynolds  Averaged  Navier-Stokes  turbulence  models  results  in  excessive  damping  of 
the  larger  scale  unsteadiness  which  is  largely  responsible  for  the  acoustic  loads  in  the  weapons 
bay.  Large  Eddy  Simulation  (LES)  resolves  all  unsteady  scales  larger  than  the  grid  related  filter 
length.  As  a  result,  not  only  does  LES  preserve  the  large  scale  acoustic  structures,  it  also  simulates 
the  largest  turbulent  scales  and  can  capture  the  interaction  between  these  scales  and  the  acoustic 
structures.  In  a  previous,  closely  related  AFOSR  program  “Large  Eddy  Simulation  for  Analysis  of 
Transonic  Cavity  Flows^,”  compressible  LES  methods  were  developed  for  application  to  cavity 
flows.  This  technology  was  implemented  in  the  LM  Aeronautics  “Falcon”  structured  grid  flow 
solver  and  applied  to  simple  cavity  flow  problems.'  Good  agreement  between  simulated  and 
experimental  acoustic  loads  was  obtained,  although  the  results  were  shown  to  be  somewhat  sensi¬ 
tive  to  grid  resolution. 

Store  release  poses  several  computational  challenges.  The  simulation  of  store  release 
requires  the  ability  to  modify  the  computational  mesh  to  capture  the  moving  store.  At  each  com¬ 
putational  iteration,  aerodynamic  loads  must  be  determined  on  the  store.  These  loads  are  com¬ 
bined  with  the  store  translation  rate,  rotation  rate,  mass  and  moment  of  inertia  in  a  six  degree  of 
freedom  dynamics  procedure  to  determine  the  displacement  of  the  store.  The  CFD  solver  must 
then  displace  the  store  computationally  by  adjusting  the  computational  mesh  to  accommodate  the 
store’s  new  position.  Effects  of  grid  motion  must  be  included  in  the  flow  solver. 

Both  structured  and  unstructured  grid  flow  solvers  have  been  used  to  simulate  store  separa¬ 
tion.  For  structured  grid  methods,  one  mesh  is  generated  around  the  store,  and  another  mesh  is 
generated  around  the  aircraft.  As  the  store  mesh  moves  with  the  store,  the  aircraft  grid  is  cut  at 
every  iteration  at  the  interface  between  the  store  grid  and  the  aircraft  grid.  At  this  interface,  mass, 
momentum  and  energy  fluxes  are  passed  between  the  two  grids.  Typically,  for  unstructured  grid 
methods,  a  single  computational  mesh  is  stretched,  distorted  and  possibly  refined  to  capture  the 
motion  of  the  store.  Typically,  if  the  store  translation  or  rotation  is  too  great,  the  unstructured  grid 
will  become  excessively  distorted,  and  generation  of  a  new  grid  is  required.  In  contrast,  the  struc¬ 
tured  grid  method  avoids  this  problem,  but  the  cutting  of  the  aircraft  grid  on  every  iteration  is 
complex  and  can  be  computationally  time  consuming. 
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1.3  Overview  of  the  Technical  Approach 

A  fighter  aircraft  with  an  internal  weapons  bay  presents  an  extreme  challenge  for  CFD  sim¬ 
ulation.  The  generation  of  a  structured  mesh  around  a  fighter  aircraft  with  an  open,  fully  loaded 
weapons  bay  would  require  many  weeks  for  an  expert  in  grid  generation.  Unstructured  mesh  gen¬ 
eration  methods  can  significantly  reduce  the  grid  generation  time  required  for  complex  configura¬ 
tions.  The  Splitflow  unstructured  grid  flow  solver  has  unique  capabilities  that  make  it  particularly 
well  suited  to  weapons  bays  analysis.  The  grid  cutting  technique  allows  complex  weapons  bay 
geometries  to  be  handled  with  relative  ease.  Starting  from  an  appropriately  processed  aircraft 
geometry,  the  Lockheed  Martin  Splitflow  unstructured  grid  flow  solver  can  be  used  to  generate  an 
initial  mesh  for  flow  solution  in  a  few  hours. 

Store  separation  simulation  with  Splitflow  is  unique,  unlike  either  structured  or  unstructured 
methods.  As  the  store  translates  through  the  flowfield,  the  Cartesian  Splitflow  mesh  is  cut  by  the 
store  geometry.  The  irregular  shaped  cells  cut  by  the  store  are  treated  as  boundary  cut  cells.  As  the 
store  translates  through  the  flowfield,  the  Splitflow  mesh  is  refined  in  areas  close  to  the  store' and 
coarsened  in  regions  which  the  store  has  passed  through.  Thus,  unlike  other  unstructured  grid 
methods,  the  grid  is  not  excessively  distorted  by  store  motion.  Splitflow  also  includes  an  alternate 
moving  body  method  where  a  body  conforming  mesh  is  generated  around  the  store,  and  this  grid 
moves  with  the  store  similar  to  the  structured  grid  store  separation  approach.  The  interfacing 
between  the  moving  store  grid  and  the  aircraft  grid  is  highly  automated  in  Spljtflow. 

Grid  resolution  is  provided  by  subdivision  of  the  Cartesian  cells.  Unlike  most  structured  and 
unstructured  flow  solvers,  individual  cells  can  be  subdivided,  while  neighboring  cells  are  left 
unchanged.  This  makes  rapid,  efficient  and  conservative  grid  refinement  possible.  This  approach 
is  extremely  well  suited  to  LES,  since  cells  can  be  refined  in  all  three  coordinate  directions  near 
walls  where  eddy  sizes  demand  high  resolution.  Away  from  walls,  cell  spacing  in  lateral  and 
streamwise  directions  is  not  governed  by  near  wall  requirements. 

1.4  Summary  of  preceding  work 

The  technical  foundation  for  the  current  program  was  developed  under  predecessor  contract 
research  and  LM  Aero  IR&D  programs.  In  the  closely  related  AFOSR  program  “Large  Eddy  Sim¬ 
ulation  for  Analysis  of  Transonic  Cavity  Flows,”  LES  capability  was  developed  for  Splitflow  and 
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implemented  in  a  hybrid  version  of  Splitflow  that  did  not  include  parallel  processing  capability. 
Between  the  completion  of  this  CRAD  program  and  the  beginning  of  the  current  research  pro¬ 
gram,  LES  modeling  capability  was  implemented  in  the  parallel  version  of  Splitflow.  In  addition 
to  parallel  processing  capability,  this  version  of  Splitflow  includes  new  grid  generation  capabili¬ 
ties  and  a  variety  of  flow  solver  enhancements. 

In  the  predecessor  program,  “Large  Eddy  Simulation  for  Analysis  of  Transonic  Cavity 
Flows."  a  compressible  form  of  the  Smagorinsky  sub  grid  scale  model  was  implemented  in  the 
hybrid  version  of  Splitflow.  In  this  version  of  Splitflow,  near  wall  boundary  layers  are  resolved 
using  triangular  prisms  generated  in  layers  by  marching  a  triangulated  surface  mesh  away  from 
the  no-slip  wall  boundaries.  The  outer  layer  of  cell  faces  from  the  prismatic  mesh  interface  with 
the  octree  based  Cartesian  Splitflow  grid.  The  Cartesian  grid  is  limited  to  cubic  cells,  and  grid 
refinement  requires  subdivision  of  a  cell  into  eight  “child”  cells.  This  approach  has  several  short¬ 
comings.  It  has  proven  difficult  to  march  the  prismatic  mesh  far  from  the  surface  for  complex 
geometries.  As  a  result,  resolution  of  the  outer  part  of  the  boundary  layer  requires  the  use  of  large 
numbers  of  cubic  cells.  In  addition,  it  is  not  possible  to  refine  the  prismatic  mesh  to  local  flow  gra¬ 
dients  because  the  location  and  density  of  the  triangular  surface  mesh  was  fixed  on  the  boundary 
surface  prior  to  running  Splitflow.  As  are  result,  automatic  refinement  capability,  a  key  feature  of 
Splitflow,  was  lost  for  the  prismatic  mesh.  The  compressible  flow  solver  was  used  to  simulate  sub¬ 
sonic.  essentially  incompressible  channel  and  backward  facing  step  flows.  Reasonable  agreement 
between  experimental  data  and  simulations  were  obtained.  Because  it  does  not  include  parallel 
processing,  the  code  took  a  large  amount  of  wall  time  to  complete  the  computations. 

Between  the  completion  of  the  LES  for  Analysis  of  Transonic  Cavity  Flows  program  and 
the  beginning  of  the  current  program,  Splitflow  LES  capability  was  significantly  improved.  The 
LES  model  was  implemented  in  a  parallel  version  of  Splitflow  with  extensive  new  capabilities. 
The  parallel  capability  allows  multiple  processors  to  be  applied  to  a  single  problem,  greatly  reduc¬ 
ing  wall  time  required  for  solutions.  The  code  can  accept  an  externally  generated,  body  conform¬ 
ing  grid  generated  by  the  “Omnigrid"  grid  generator.  Grids  generated  with  Omnigrid  can  be 
refined  in  a  manner  similar  to  the  traditional,  “cut  grid”  version  of  Splitflow.  This  version  of  Split- 
flow  can  also  generate  its  own  Cartesian  mesh.  Unlike  the  hybrid  version  of  Splitflow,  the  Carte¬ 
sian  mesh  is  not  restricted  to  cubic  cells.  Cells  can  be  split  in  any  of  the  3  Cartesian  directions. 
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creating  either  two,  four  or  eight  children  in  a  refinement  process,  resulting  in  non-cubic  cells. 
Aircraft  geometries  and  flowfield  gradients  can  be  resolved  with  far  fewer  cells  than  was  required 
with  the  oct-tree  version  of  Splitflow. 


Section  2:  Methodology 

In  order  to  present  a  coherent  picture  of  the  methodology  applied  to  the  store  separation  and 
acoustic  problems,  methods  developed  prior  to  the  current  program  as  well  as  methods  developed 
as  part  of  the  current  program  are  described  in  Section  2.  Section  2.6,  LES  Modeling  and  Imple¬ 
mentation,  describes  methods  developed  under  the  “Large  Eddy  Simulation  for  the  Analysis  of 
Transonic  Cavity  Flows”  program  and  implemented  in  the  parallel  version  of  Splitflow  under  an 
LM  Aero  IR&D  program.  Section  2.1,  Parallel  Splitflow  and  Section  2.2,  Splitflow  Grid  Schemes 
describe  methods  developed  and  implemented  under  IR&D  prior  to  the  beginning  of  the  current 
program.  The  balance  of  Section  2  describes  methods  developed  as  part  of  the  current  program. 

2.1  Parallel  Splitflow 

Parallel  processing  is  achieved  using  the  Parallel  Virtual  Machine  (PVM)  message  passing 
protocol.  One  process  of  each  parallel  job  is  designated  as  the  “mother”  process,  where  grid  gen¬ 
eration,  post  processing  and  the  other  parallel  processes  are  controlled.  The  grid  encompassing 
the  entire  flow  domain  is  decomposed  into  a  number  of  subdomains,  each  of  which  allocated  to  a 
proc.'.ss.  The  subdomains  are  defined  automatically,  and  the  number  of  cells  and  number  of 
boundary  triangular  facets  can  be  balanced  in  order  to  obtain  efficient  use  of  multiple  processors. 

A  pointwise  implicit  solution  algorithm  is  usually  used  for  parallel  processing  in  Splitflow. 
The  pointwise  scheme  iteration  history  is  independent  of  the  number  of  processors  applied.  For 
the  pointwise  scheme,  multiple  subiterations  are  performed  at  each  iteration  level  to  obtain  a  time 
accurate  solution.  The  basic  numerical  form  of  the  implicit  equation  is: 
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where  the  superscript  n  indicates  the  previous  time  step.  The  superscript  s  designates  the  most 
updated  estimate  of  time  level  n+J.  When  the  solution  at  sub-iteration  s+1  is  sufficiently  equal  to 
that  at  sub-iteration  s,  then  the  sub-iterations  are  complete,  and  the  solution  is  advanced  to  n+1. 

The  flux  Jacobians  include  the  inviscid  Jacobians  consistent  with  Roe’s  scheme,  computed 
using  first  order  data.  The  algorithm  requires  a  block  inversion  of  a  5  x  5  matrix  for  each  cell.  The 
inverted  matrix  is  computed  during  the  first  sub-iteration  and  stored  for  use  in  subsequent  sub¬ 
iterations.  Although  the  Jacobians  are  frozen,  the  fluxes  must  be  updated  with  each  sub-iteration. 
Upon  convergence  of  the  subiterations,  the  linearization  errors  are  eliminated  because  the  inter¬ 
face  fluxes  are  evaluated  directly  from  the  most  up-to-date  solution  at  the  centroids.  Typically,  10 
to  1 5  sub-iterations  are  adequate  to  sufficiently  converge  the  implicit  equations  for  a  time  depen¬ 
dent  solution.  More  sub-iterations  are  generally  required  for  as  the  time  step  is  increased.  A  more 
detailed  description  of  Splitflow  solution  algorithms  is  contained  in  Domel  and  Karman.^ 

2.2  Splitflow  Grid  Schemes 

Cut  Grid  Method:  In  the  cut  grid  scheme,  cell  faces  are  oriented  with  respect  to  Cartesian 
axes.  Boundary  surfaces,  including  no-slip  walls,  are  represented  by  triangulated  surface  meshes. 
Where  the  triangulated  surface  mesh  intersects  Cartesian  cells,  the  Cartesian  cells  are  cut  into 
irregular  shapes.  The  cut  cell  surfaces  are  cut  into  triangular  boundary  facets  as  shown  in  Figure- 
Equation  2.1.  Boundary  conditions  are  applied  at  the  boundary  facets,  and  appropriate  fluxes  are 
evaluated  at  each  boundary  cut  cell.  Because  of  this  approach,  the  Splitflow  simulation  can  com¬ 
mence  once  a  triangulated  surface  mesh  is  defined. 


Figure  2.1  Splitflow  automated  cut  grid  method  allows  rapid  problem  set-up 

Grid  refinement  is  achieved  by  cell  sub-division.  Cells  can  be  divided  in  any  of  three  coordi¬ 
nate  directions.  Figure  2.15  shows  the  cell  division  process.  Cells  are  divided  in  an  initial  grid 
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generation  process  to  meet  geometry  resolution  constraints.  Typically,  10-20  levels  of  cell  divi¬ 
sion  are  present  in  an  initial  grid.  Once  the  flow  solver  is  applied  to  the  initial  grid,  addition  cells 
are  refined  in  order  to  resolve  flowfield  gradients. 


2  children  4  children 


2  children  4  children 


2  children  4  children 

Figure  2.2  Splitflow  grid  cell  refinement  options 

“Omnigrid”  Body  Fitted  Grid  Method:  Omnigrid  is  used  to  develop  body  conforming  grids 
with  sufficient  smoothness  and  near  wall  resolution  for  efficient  viscous  flow  simulations.  The 
omnigrid  process  begins  with  generation  of  a  conventional  oct-tree  based  Cartesian  Splitflow 
mesh  around  a  geometry  represented  by  a  triangulated  surface  facets.  Next,  cells  intersected  by  or 
in  close  proximity  to  the  geometry  are  removed.  New  grid  lines  are  projected  from  the  nodes  on 
the  exposed  faces  of  the  remaining  cells  to  the  surface.  This  produces  a  body  fitted  grid.  High  res¬ 
olution  grids  for  viscous  analysis  are  generated  by  refining  the  grid,  particularly  in  the  layer  of 
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cells  between  the  body  surface  and  the  Cartesian  volume  mesh.  The  resulting  mesh  can  be  refined 
to  resolve  other  flow  field  gradients  using  user  selected  adaptation  functions. 

Resolution  of  sharp  features  in  the  geometry  can  be  critical  for  accurate  simulation.  Omni¬ 
grid  can  capture  sharp  features  in  the  geometry  by  “snapping”  node  points  to  edge  locations.  To 
preserve  sharp  edges  in  the  geometry,  quadrilateral  faces  on  the  surface  of  the  geometry  in  the 
vicinity  of  the  sharp  edge  are  identified.  Using  these  faces,  the  mesh  nodes  closest  to  the  sharp 
edge  are  labeled  and  snapped  to  sharp  edges.  This  approach  allows  ch'e  computational  mesh  to 
capture  geometric  details  without  excessive  numbers  of  nodes. 

The  omnigrid  generation  process  in  general,  and  the  snapping  process  in  particular,  can  cre¬ 
ate  regions  of  tangled  mesh  with  grid  crossing.  To  correct  these  regions  and  improve  the  quality  of 
the  grid,  an  elliptic  smoothing  scheme  is  used  for  interior  points  with  a  spring  analogy  scheme  for 
boundary  nodes.  The  current  version  of  the  elliptic  smoothing  solves  the  Laplace  equations.  An 
example  of  the  smoothing  scheme  is  shown  in  Figure  2.3  for  the  cross  section  of  a  chined  fore¬ 
body. 


Figure  2.3  Example  of  a  grid  around  a  chined  nozzle  generated  with  Omnigrid 
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2.3  Moving  Grid  Methods 

The  addition  of  rigid  body  motion  complicates  the  time-accurate  CFD  simulation  by  intro¬ 
ducing  grid  speed  terms  and  changing  cell  volumes.  Investigators  generally  approach  this  com¬ 
plexity  in  one  (or  some  combination)  of  three  ways,  each  having  advantages  and  disadvantages. 
The  first  method  involves  a  single  deforming  grid  which  preserves  the  identity  of  each  grid  cell. 
Each  cell  is  allowed  to  stretch  and  deform  in  order  to  fill  the  computational  domain  around  the 
moving  body.  This  works  well  for  small  displacements,  but  large  displacements  often  result  in 
highly  skewed,  compressed  or  stretched  grid  cells  which  necessitate  occasional  regridding  and 
redistributing  cells. 

The  second  method,  overset  grids,  allows  rigid  grids  to  be  defined,  thereby  maintaining  the 
initial  grid  resolution.  An  overset  grid  moves  with  a  body  over  a  stationary  background  grid.  The 
background  grid  has  a  “hole”  which  is  larger  than  the  moving  body,  but  smaller  than  the  moving 
grid  surrounding  that  body.  The  boundaries  of  the  “hole”  interface  with  the  solution  in  the  moving 
grid.  Similarly,  the  outer  boundaries  of  the  moving  grid  interface  with  the  background  grid.  Cells 
in  the  background  grid  are  turned  off  and  on  as  needed  in  order  to  allow  the  hole  to  follow  the 
moving  body.  The  extent  of  the  overlap  between  the  background  grid  and  the  moving  grid  must 
contain  the  entire  motion  of  the  body  from  one  time  step  to  the  next,  thereby  ensuring  that  each 
point  in  the  computational  region  is  defined  in  at  least  one  of  the  grids.  Grids  with  large  overlaps 
allow  large  displacements  between  time  steps.  However,  care  must  be  taken  as  the  moving  body 
nears  the  boundary  of  the  background  grid  because  the  moving  grid  may  extend  into  undefined 
territory  outside  the  computational  domain. 

The  third  approach  is  to  redefine  boundary  cells  every  iteration  as  the  body  moves.  This 
approach  allows  for  the  majority  of  the  cells  to  remain  fixed  from  one  time  step  to  the  next,  with 
only  boundary  cells  changing  with  time.  This  approach  allows  for  very  extreme  motion  of  the 
body  within  the  computational  domain,  but  the  appearance  and  disappearance  of  cells  from  within 
the  computational  domain  must  be  addressed.  This  third  method  was  adopted  for  Splitfiow 
because  of  its  automatic  grid  generation  capability,  and  its  ability  to  modify  a  pre-existing  grid  to 
contain  a  new  geometry  via  the  cell  cutting  process. 
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Splitflow  models  the  motion  of  the  geometry  by  updating  the  position  of  the  store  at  every 
iteration,  and  recutting  the  Cartesian  grid  in  order  to  accommodate  the  updated  geometry. 
Although  the  underlying  Cartesian  mesh  does  not  move  with  time,  some  Cartesian  cells  may 
undergo  refinement  in  order  to  resolve  the  new  geometry  position.  After  recutting  the  Cartesian 
cells  to  fit  around  the  new  geometry  position,  Splitflow  assigns  grid  speeds  to  the  appropriate 
boundary  facets.  The  grid  speeds  of  these  boundary  facets  determine  the  rate  of  change  of  the  vol¬ 
umes  of  the  boundary  cells  adjacent  to  the  moving  body.  Only  cells  adjacent  to  a  moving  bound¬ 
ary  have  time-variant  volumes. 

This  method  of  recutting  the  grid  around  a  moving  body  causes  new  cells  to  appear  inside 
the  computational  domain.  Similarly,  cells  may  disappear  from  the  computational  domain  due  to 
the  motion  of  the  body.  The  proper  way  to  strictly  enforce  conservation  in  such  cells  is  to  merge 
each  one  with  neighboring  cells  such  that  each  aggregate  merged  cell  simply  changes  volume, 
rather  than  appearing  or  disappearing  completely.  However,  the  algorithm  for  determining  such 
cell  combinations  may  get  quite  complicated  if  the  moving  boundary  is  allowed  to  traverse  several 
cells  in  one  iteration.  For  this  program,  cells  are  allowed  to  appear  and  disappear.  While  this 
approach  is  not  strictly  conservative,  it  provides  adequate  accuracy  and  reduces  computational 
solution  times.  Cells  are  allowed  to  disappear  without  attention.  However,  attention  must  be 
focu.sed  upon  the  initialization  of  the  solution  in  new  cells  which  appear  in  the  computational 
domain. 

A  method  was  developed  for  Splitflow  which  approximates  the  solution  from  an  overset 
grid  moving  with  the  body.  Because  an  actual  overset  grid  does  not  exist,  Splitflow  estimates  var¬ 
ious  parts  of  the  solution  as  needed.  Splitflow  calculates  the  position  of  each  uninitialized  cell  rel¬ 
ative  to  the  moving  body.  Then  Splitflow  interpolates  from  the  appropriate  position  in  the  solution 
at  the  previous  geometry  location,  taking  into  account  the  translation  and  rotation  since  the  previ¬ 
ous  lime  level.  Figure  2.4  illustrates  the  method. 

Although  the  described  method  provides  good  accuracy  for  the  simulation,  the  computa¬ 
tional  time  required  for  searching  and  interpolating  warranted  a  compromise  for  faster  initializa¬ 
tion  from  ancestor  cell  information.  Splitflow  accomplishes  this  by  recognizing  that  an 
“appearing”  cell  is  actually  a  pre-existing  cell  which  switches  its  state  from  “off’  to  “on.”  There¬ 
fore.  Splitflow  determines  a  solution  for  each  “off’  cell  before  it  switches  state.  Before  the  body 
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position  is  updated  for  an  iteration,  Splitflow  assigns  each  parent  cell  a  solution  based  upon  the 
volume  weighted  average  of  its  computational  children.  Splitflow  systematically  works  its  way  up 
to  the  root  cell,  assigning  a  solution  to  every  parent  in  the  computational  domain.  Then  Splitflow 
reverses  direction  and  assigns  solutions  to  each  cell  outside  the  computational  domain  based  upon 
its  nearest  computational  ancestor.  Therefore,  all  cells  outside  the  computational  domain  may  be 
initialized  prior  to  moving  the  body  for  any  iteration. 


Figure  2.4  Method  fOi_  initializing  cells  which  appear  due  to  geometry  motion 


One  other  option  included  with  the  moving  body  computations  in  Splitflow  is  the  capability 
to  freeze  the  body  position  (and  grid  speeds  and  volume  rates)  for  a  few  iterations,  and  then  com¬ 
pensate  when  the  body  position  is  actually  updated.  This  reduces  frequency  of  body  motion  and 
grid  cutting,  which  can  significantly  reduce  solution  times.  Some  sacrifice  in  solution  accuracy  is 
possible  because  the  body  position  is  updated  less  frequently,  resulting  in  greater  error  in  store 
positioning.  However,  the  robustness  corresponds  to  that  of  the  smaller  time  step  between  itera¬ 
tions. 
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2.4  Interface  Grid  Methods 

Flow  field  data  transfer  between  domains  was  accomplished  through  an  interfacing  system 
developed  at  LM  Aero.  This  flexible  interfacing  method  uses  a  object-oriented  file  structure  to 
pass  flow  field  data  and,  in  some  cases,  flux  data  between  otherwise  unconnected  domains.  It  can 
also  be  used  to  pass  information  between  separate  Splitflow  instances,  or  processes,  and  allow  for 
simultaneous  solving  on  mixed  grid  types.  A  description  of  the  interfacing  method  and  how  it  was 
used  for  this  program  follows. 

Method 

Two  types  of  interfaces  have  been  developed.  The  type  1  interface  is  used  to  communicate 
between  two  Splitflow  cut-type  grids.  Interface  type  2  can  communicate  between  a  Splitflow  cut 
grid  and  a  body-fitted  grid  (Figure  2.5).  The  body-fitted  grid  could  theoretically  be  any  structured 
or  unstructured  grid.  Both  interface  types  share  flow  field  data  across  the  interface.  Flow  field  data 
includes  density,  three  momentum  components,  and  energy  and  will  hereafter  be  referred  to  as 
“Q”  data.  Flow  field  data  are  obtained  from  the  solutions  on  each  side  of  the  interface.  The  Roe 
scheme  is  used  to  evaluate  the  boundary  fluxes  that  are  applied  to  both  solutions.  The  resulting 
data  transfer  conserves  mass,  momentum,  and  energy. 

Identical  geometry  triangulation  is  required  for  both  domains  so  that  a  common  set  of  areas 
and  orderings  can  be  assumed.  The  resolution  or  topology  of  this  triangulation  is  not  restricted  for 

Type  1  Interface  Type  2  Interface 


Figure  2.5  Two  interface  types  are  available  in  Splitflow 
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the  type  1  interface  but  resolution  should  be  on  the  same  order  as  the  associated  boundary  cells. 
For  a  type  2  interface,  the  triangulation  is  generated  by  connecting  opposite  comers  of  the  body- 
fitted  grid  quadrilateral  elements  and  using  that  on  the  cut  grid  side.  Interfacing  with  a  tetrahedral 
grid  could  be  done  without  modifying  the  surface  mesh. 

The  sequence  for  a  type  1  interface  on  Domain  1  (Figure  2.6)  is  as  follows: 

1 .  Domain  1  collects  local  boundary  solution  to  the  geometry  facets. 

2.  Local  geometry  facet  solutions  are  exported  to  Domain  2. 

3.3)  Geometry  facet  data  is  imported  from  Domain  2. 

4.  The  Roe  scheme  is  applied  to  each  geometry  facet  using  the  solutions  from  both  sides  of  the 
interface. 

5.  The  resulting  fluxes  are  applied  to  both  domains. 

The  same  sequence  takes  place  on  Domain  2  so  that  the  same  Q  and  flux  information  is  applied  to 
both  sides  of  the  boundary. 


Domain  1 
Cut  Grid 


‘local  (export)^ 


i-Q 


Roe 


import 


t 

Rux 


‘local  (export) 


TV 

Q 


import^ 


Roe 


►1*^ - Flux 


Domain  2 
Cut  Grid 


Figure  2.6  Type  1  interface  computation  and  transfer  sequence 

The  sequence  for  a  type  2  interface  (Figure  2.7)  is  as  follows: 

1 .  Domains  1  and  2  collect  local  boundary  solution  to  the  geometry  facets. 

2.  Domain  2  exports  local  boundary  solution  to  Domain  1. 

3.  Domain  1  averages  the  solutions  (Qbound)  applies  Roe  scheme  to  obtain  fluxes. 

4.  Boundary  Q  data  and  fluxes  are  applied  in  Domain  1 
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5.  Boundary  Q  and  flux  data  are  transferred  to  Domain  2  (Qbound,  Flux) 

6.  Boundary'  Q  and  flux  data  are  applied  in  Domain  2 


Domain  1 
Cut  Grid 


■local 


y  ^import 


—  L 

i 


■bound 


Roe 


Flux 


■local{export) 


Domain  2 
Body-Fitted  Grid 


Figure  2.7  Type  2  interface  computation  and  transfer  sequence 

Data  can  be  passed  in  one  direction  as  well.  This  amounts  to  skipping  step  2  and  using  Qj 
data  on  both  sides  of  the  Roe  scheme.  A  similar  method  can  also  be  used  to  get  data  from  Domain 
1  at  some  point  in  time  to  be  used  as  a  non-uniform  fixed  boundary  input  for  Domain  2.  The  “out¬ 
flowing”  domain  can  also  apply  another  boundary  condition  in  conjunction  with  the  interface.  For 
example,  an  extrapolation  boundary  condition  could  be  applied  and  those  Q  data  would  also  be 
transmitted,  through  the  interface,  to  another  domain. 


Synchronization  can  be  used  to  keep  both  solutions  advancing  in  time  at  the  same  rate.  Thir 
is  done  on  the  Iteration  level  for  local  time  stepping  (time  warped)  and  on  the  sub-iteration  level 
for  global  time  stepping  (time  accurate).  Rotations  can  also  be  applied  to  the  momentum  compo¬ 
nents  at  the  interface  to  account  for  arbitrary  relative  grid  orientations.  This  might  be  used  for  a 
case  where  one  grid  is  moving  (with  rotations)  through  another. 


The  communication  path  for  these  interfaces  is  a  file  that  is  accessible  to  both  domains.  A 
portable,  self-de.scribing,  object-oriented  file  format  is  used  which  was  developed  by  LM  Aero. 
Such  an  object-oriented  approach  allows  for  a  theoretically  infinite  number  of  interfaces,  although 
file  input/output  (I/O)  times,  etc.  impose  a  practical  limit.  On  start-up,  Splitflow  calculates  inter¬ 
face  sizes  and  allocates  the  necessary  memory.  The  first  time  communication  is  attempted,  exist¬ 
ence  of  the  interface  file  is  checked.  If  no  file  exists,  one  is  created  with  an  (I/O)  blocking  object 
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and  interface  objects  are  added  as  needed.  The  I/O  blocking  object  includes  a  flag  which  gets 
turned  on  whenever  a  process  is  writing  data  into  the  file.  This  allows  simultaneous  reading  but 
not  simultaneous  writing  or  reading  while  another  process  is  writing,  since  the  latter  two  situa¬ 
tions  could  result  in  corrupted  data.  In  the  future  a  message  passing  protocol  could  be  used  in 

place  of  the  interface  file  method. 

> 

Simple  Cavity  Application 

For  the  simple  cavity  calculation,  two  interfaces  were  used.  One  interface  was  used  at  the 
outflow  of  the  channel  domain  to  “copy”  data  to  a  transition  domain.  The  other  linked  the  transi¬ 
tion  domain  to  the  cavity  domain  (Figure  2.8).  Since  the  channel  and  transition  domains  resided  in 
the  same  grid  and  were  run  as  one  instance  of  the  Splitflow  solver,  no  synchronization  was  neces¬ 
sary.  The  synchronization  feature  was  used  for  the  type  2  interface  to  keep  the  two  solutions  at  the 
same  time  level. 


Type  2  Interface 


F-22  Cavity  Application 

To  produce  an  appropriate  incoming  boundary  layer  for  the  F-22  main  weapons  bay,  a 
body-fitted  grid  was  generated  for  the  part  of  the  F-22  forebody.  This  forebody  grid  communi¬ 
cated  with  a  cut  grid  defining  the  rest  of  the  aircraft  through  an  irregularly-shaped  interface 
boundary.  Figure  2.9  shows  the  aircraft  surface  boundaries  included  in  each  of  the  domains  -  gray 
in  the  cut  grid  and  cyan  in  the  body-fitted  grid.  In  Figure  2.10  the  interface  boundary  is  shown.  A 
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Gray  -  Cut  grid  Cyan  -  body-fitted  grid 


Figure  2.9  Aircraft  surfaces  solved  in  each  domain 


Front  View  Side  View  Bottom  View 


Figure  2.10  The  forebody  interface 

two-direction,  type  2  interface  was  applied  on  this  boundary  to  pass  Q  and  flux  data  between  the 
grids  (Figure  2.1 1  and  Figure  2.12).  The  two  grids  were  run  as  separate  instances  of  the  Splitflow 
solver  and  synchronized  on  the  sub-iteration  level. 


Figure  2.11  Fuselage  station  cuts  through  both  grids 
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Figure  2.12  Butt  line  cut  through  both  grids  with  close-up  side  view 


2.5  Six  Degree  of  Freedom  Methods 

Following  the  addition  of  time-accurate  body  motion  in  Splitflow,  six  degree-of-freedom 
(6DOF)  capability  was  implemented.  Constant-rate  motion  was  first  introduced  and  tested  to 
show  that  the  grid  speed  terms  were  properly  implemented.  Having  accomplished  this,  body 
motion  capability  was  expanded  to  include  arbitrary  motion.  Two  types  of  arbitrary  motion  are 
allowed  for:  1 )  prescribed  and  2)  6DOF.  A  restraint  module  is  also  available  in  the  6DOF  routines 
which  can  be  used  to  reduce  degrees  of  freedom  as  desired. 

Assumptions  and  Definitions 

The  6DOF  routines  in  Splitflow  use  the  Euler  equations  of  motion  as  described  in  Dynamics 
of  Flight:  Stability  and  Control,  second  edition,  by  Bernard  Etkin,  P.  88.'^  From  inertial  properties, 
forces  and  moments,  and  initial  linear  and  rotational  velocities,  linear  and  rotational  accelerations 
can  be  computed.  An  integration  in  time  produces  linear  and  rotational  velocities.  Another  inte¬ 
gration  in  time  produces  position  and  orientation. 

Store  Separation 

The  6DOF  routines  implemented  in  Splitflow  were  originally  developed  for  store  separation 
simulation.  As  a  result,  many  of  the  conventions  used  are  those  which  are  common  to  store  sepa¬ 
ration  analyses.  Axis  system  conventions  (described  in  a  following  section)  are  one  of  the  more 
imponant  outcomes  of  this  approach  in  that  it  becomes  important  for  the  user  to  understand  them 
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to  properly  set  up  a  simulation.  The  inclusion  of  aircraft  maneuver  (restricted  to  the  pitch  plane, 
i.e.  pull-up  and  push  over)  and  ejector  modeling  are  other  features  useful  for  store  separation  that 
may  not  be  of  value  in  more  generalized  use. 

Axis  System  Convention 

Aircraft  geometry  is  typically  defined  in  a  fuselage  station,  bur.  line,  water  line  system 
(upper  left  in  Figure  2.13).  The  convention  commonly  used  for  store  separation  analysis  is  called 
a  body  axis  system.  This  system,  along  with  the  corresponding  force  and  moment  coefficient  def¬ 
initions,  is  shown  on  the  right  in  Figure  2.13.  Transforming  from  a  fuselage  station,  butt  line, 
water  line  system  to  a  body  system  is  simply  a  matter  of  rotating  180  degrees  about  the  Y  axis. 
Origins  of  either  of  these  systems  are  typically  located  forward  and  below  the  nose  of  the  aircraft 
on  centerline  or  at  the  aircraft  center  of  gravity  (CG).  A  pylon  axis  system  (lower  left  in 
Figure  2.13)  is  defined  by  a  store  body  axis  system  at  the  carriage  position.  Most  trajectory  data 
are  provided  in  a  pylon  system. 


Figure  2.13  6DOF  axis  system  conventions 


Integration  Scheme 

A  trapezoidal  scheme  is  used  to  integrate  accelerations  from  the  Euler  equations  of  motion 
to  obtain  velocities.  Another  integration  is  performed  to  produce  position  and  orientation.  Higher 
order  integration  schemes  have  been  investigated  in  the  past  as  part  of  an  effort  to  reduce  trajec- 
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tory  error.  The  result  of  the  study  was  that,  since  these  routines  are  not  computationally  intensive, 
simply  reducing  the  integration  time  step  was  the  easiest  and  most  effective  means  of  improving 
accuracy.  As  implemented  in  Splitflow,  the  integration  time  step  is  independent  of  the  solution 
time  step  so  an  integration  time  step  one  to  two  orders  of  magnitude  smaller  than  the  solution  time 
step  is  recommended  (something  on  the  order  of  1.0x10’^  or  1.0x10-7  seconds  per  integration 
step) 

Models  and  Initial  Conditions 

Aircraft  maneuver  modeling  is  limited  to  the  pitch  plane  (pull-up  and  push-over).  Angular 
velocity  of  the  aircraft  is  computed  from: 


Q  =  (2) 

where  Q  is  angular  velocity  in  the  Y-direction,  n  is  the  load  factor,  g  is  acceleration  due  to 
gravity,  and  V  is  aircraft  velocity  tangent  to  the  maneuver  radius.  A  load  factor  of  one  corresponds 
to  straight  and  level  flight  (no  maneuver),  a  load  factor  greater  than  one  corresponds  to  a  pull-up 
maneuver,  and  a  load  factor  less  than  one  corresponds  to  a  push-over  maneuver. 

Two  store  ejectors  can  be  modeled  by  supplying  ejector  positions,  end  of  stroke  (EOS) 
length,  and  a  force-time  history  of  the  ejector.  Typical  ejectors  produce  a  somewhat  “sawtooth” 
force-time  history  over  about  60  milliseconds  with  enough  force  to  accelerate  the  store  to  ~30  feet 
per  second  at  EOS.  Often  the  forward  ejector  force  is  higher  to  produce  a  nose-down  angular 
velocity.  An  alternative  to  ejector  modeling  is  setting  initial  conditions,  such  as  EOS  linear  and 
angular  velocities. 

6DOF  Code  Flow  Chart 

SOLVE 

Call  MOVE_BODIES 

•Solve  flow  field  for  current  iteration 
•Sum  aerodynamic  loads  on  the  body 

•Append  current  body  position,  orientation,  and  aerodynamic  loads  to  trajectoiy  history  file 
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END 

MOVE_BODIES 

•Collect  boundaries  acting  as  one  arbitrarily  moving  body 
•Sum  aerodynamic  forces  and  moments  on  the  body 

Call  ARB_M0T10N 

•Initialize  newly  uncovered  cells 
•Move  geometry 
•Recut  grid 

•Update  moment  center 

RETURN 

ARB_MOTION 

•Read  input  file 

•Attempt  to  open  trajectory  history  file  and  read  trajectory  data 
•Scale  body  loads 

•Transform  body  loads  into  the  inertial  axis  system 
•Transform  body  loads  into  the  body  axis  system 
•Establish  end  time 

Call  SIXDOFJNTEG 

•Scale  new  body  position,  orientation,  and  rotational  rates 

RETURN 

SIXDOF_INTEG 

•Set  up  aircraft  motion  terms  (maneuver  in  pitch  plane) 

•Set  up  store  motion  terms  as  a  result  of  aircraft  maneuver 
•Start  integration  loop 

••Interpolate  store  loads  in  Z 

••Interpolate  other  inputs  in  time  (ejector,  etc.) 
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••Formulate  acceleration  terms 
••Integrate  to  get  velocities 
••Integrate  to  get  position  and  orientation 
••Transform  results  into  inertial  axis  system 
••Record  results  to  output  files 
••Check  for  stop  condition 

RETURN 


2.6  LES  Modeling  and  Implementation 

LES  Model:  The  model  is  a  compressible  form  of  the  Smagorinsky  model.^  By  applying  a  “filter 
function”  to  the  Navier-Stokes  equations,  the  LES  form  of  the  flow  equations  can  be  derived.  The 
filter  is  a  grid  dependent  integral  of  the  flow  variables  that  eliminates  variations  in  the  flow  quanti¬ 
ties  that  are  smaller  than  can  be  resolved  by  the  computational  mesh.  A  filtered  variable,  /,  is  the 
spatial  average  using  a  filter  function,  G,  whose  integral  at  any  x  is  unity.  Favre  averaging  is  uti¬ 
lized  inside  the  filtering  operation  to  account  for  density  variations,  which  gives  the  mass  aver¬ 
aged,  filtered  variable  / . 

fix)  =  pjp/(x)G(x-y)Jy  (3) 


For  finite  volume  applications,  the  “top  hat”  filter  is  typically  used,  where: 


GiXi) 


fl/(Ax,.) 

I  0 


ifi\x-\  <  Ax/2 
ifi\xi\  >  Ax/2 


(4) 


and  Xj  is  the  component  of  the  grid  spacing  in  computational  space,  and  Ax,-  is  the  cell  spacing  in 
the  i-direction.  Applying  the  filter  function  to  the  continuity,  momentum  and  energy  equations 
gives: 


dt  dx: 


ipu.)  =  0 


(5) 
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3pw,  3  ^  5  p 


(6) 


The  subgrid  scale  momentum  stress  is  o-j  =  pu-Uj  -  pu-u j 


/  a  _  ,  ^6,1 


^dx.  ax,.^ 


=  0 


(7) 


The  subgrid  scale  diffusion  terms  are: 


9/  = 

\|/.  =  ^ip^^^^k-pUiUkUi,) 

0.  =  (pe  +  P)u--{pe  +  P)Ui 


(8) 


The  subgrid  scale  stress  term  is  modeled  as  in  Smagorinsky  with  a  simple  Bousinesq  approxima¬ 
tion. 


where 


5  _  1 

-  2 


f duj  du  ^ 


JJ 


(9) 


(10) 


The  eddy  viscosity  is: 


li,  =  p/'|5|  |5|  =  (11) 

Away  from  no  slip  walls  the  mixing  length,  /,  is  simply  a  function  of  the  grid  scale,  A,  and  the 
Smagorinsky  constant,  Cj,  /  =  C^A.  The  constant  Q  is  set  to  0.1  for  this  work.  A  simple  model 
is  also  used  for  the  subgrid  scale  heat  flux. 
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Pr^  ^.v• 


(12) 


Where  Cp  is  the  specific  heat  at  constant  pressure  and  we  set  the  turbulent  Prandtl  number,  Pr,  = 
0.9.  The  final  model  term  for  the  subgrid  scale  viscous  work  in  the  energy  equation  is: 

vi/,  +  (p.  =  UjG^j  (13) 

Given  the  value  of  the  filter  length  scale,  /,  the  continuity,  momentum  and  energy  equations.  Equa¬ 
tions  5,  6  and  7  are  closed  with  the  model  terms.  Equations  9,  12  and  13. 

Filter  Length  for  Splitflow  Grid  LES  modeling  requires  gradual  changes  in  the  filter  length 
throughout  the  flowfield  in  order  to  be  able  to  simulate  large  turbulent  structures.  In  LES,  the  filter 
length  is  generally  determined  as  a  function  of  the  local  cell  volume: 

A  =  (cell  volume)^ (14) 

Unlike  a  structured  grid,  the  a  Splitflow  mesh  can  have  neighboring  cells  with  large  disparities  in 
cell  size  due  to  cell  cutting  by  geometry  surfaces  or  grid  refinement.  The  volumes  of  two  uncut 
cells  sharing  a  common  face  will  never  differ  by  more  than  a  factor  of  eight.  Rapid  changes  in  fil¬ 
ter  length  are  inconsistent  with  sub  grid  scale  approximation.  The  small  scales  of  turbulence  will 
not  be  produced  or  simulated  well  in  a  limited  region  of  the  flow  fiel4  with  a  fine  grid  and  small 
filter  lengths  surrounded  by  coarser  grids  and  larger  filter  lengths.  Cuiy  larger  scale  structures  will 
be  convected  through  the  region  of  reduced  scales;  there  is  insufficient  time  for  the  smaller  scale 
turbulence  to  develop.  If  the  filter  length  is  not  smoothed,  the  turbulent  viscosity  will  be  much 
smaller  in  the  fine  grid  region.  This  can  cause  “kinks”  in  mean  velocity  profiles  obtained  from 
averaged  LES  solutions.  Conversely,  if  the  filter  length  in  local  regions  of  high  grid  refinement  is 
much  larger  than  the  local  grid  scale,  but  roughly  the  same  magnitude  as  adjacent  regions  with 
coarser  grids,  the  numerical  accuracy  of  the  solution  is  improved  in  the  region  of  highly  refined 
grid,  and  the  turbulent  viscosity  varies  smoothly  throughout  this  region.  For  this  reason  a  smooth¬ 
ing  scheme  that  increases  the  filter  length  smoothly  across  regions  with  rapid  changes  in  cell  vol¬ 
ume  was  devised. 
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For  each  sweep  of  the  smoothing  function,  the  maxima  of  the  filter  length  of  each  cell  and 
the  filter  length  of  its  neighboring  cells  are  averaged.  For  each  smoothing  sweep,  the  filter  length 
for  cell  at  c,  is: 


N 

n  =  1 


where,  N  is  the  number  of  neighbors.  For  all  of  the  results  presented  here,  three  smoothing  passes 
were  employed.  With  this  approach  a  reasonably  smooth  filter  length  distribution  was  obtained. 

The  Smagorinsky  constant  is  reduced  near  no  slip  wall  boundaries  using  van  Driest  damp¬ 
ing.  In  addition,  the  filter  length  is  not  allowed  to  exceed  the  product  of  the  distance  from  the  wall 
and  von  Karman  constant,  K.  Most  LES  simulations  with  Smagorinsky’s  model  do  not  make  this 
restriction,  resulting  in  poor  near  wall  results  when  large  spanwise  and  streamwise  grid  spacings 
cause  the  filter  length  to  exceed  Ky  outside  the  viscous  sublayer.  The  length  scale  is: 


/  =  min(C^A,  kDv) 


1  -  exp 


(16) 


The  constant  is  set  to  25.0.  The  friction  velocity,  is  calculated  from  a  combination  of  local 
conditions  and  wall  conditions. 


= 


(17) 


The  nearest  boundary  facet  is  identified  for  each  cell.  The  wall  shear  stress  associated  with  this 
boundary  facet  is  applied  in  Equation  17.  The  length  scale,  - 1  is  lagged  from  the  previous  itera¬ 
tion.  The  second  term  in  Equation  17  prevents  the  damping  term  from  going  to  zero  away  from 
the  wall  when  the  instantaneous  wall  shear  stress  is  zero.  This  is  reasonable  since  we  expect  the 
near  wall  damping  to  reduce  turbulent  viscosity  when  the  size  of  sub  grid  scale  turbulent  struc¬ 
tures  is  greater  than  the  distance  from  the  wall.  When  the  wall  shear  stress  is  zero,  velocity  gradi¬ 
ents  and  subgrid  scale  stresses  away  from  the  wall  are  still  significant  and  should  be  damped. 
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Wall  function  boundary  condition:  The  wall  function  boundary  condition  for  the  LES  model 
implemented  in  Splitflow  solves  for  a  velocity  profile  that  matches  the  velocity  at  the  cell  center 
one  point  away  from  the  wall  and  the  shear  stress  at  the  face  of  the  grid  cell  one  grid  cell  height 
away  from  the  wall.  This  profile  is  used  to  obtain  the  wall  shear  stress,  which  is  used  in  the  vis¬ 
cous  flux  routine.  The  wall  function  greatly  relaxes  near  wall  grid  resolution  requirements.  The 
LES  wall  function  is  more  sophisticated  than  the  law  of  the  wall  because  the  velocity  profile  con¬ 
forms  to  the  instantaneous  shear  stress  profile  between  the  wall  and  the  first  grid  point  off  of  the 
wall,  while  the  law  of  the  wall  assumes  constant  shear  stress. 
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Figure  2.14  Notation  for  wall  function  grid  locations. 

A  diagram  of  the  wall  function  locations  is  shown  in  Figure  2.14.  The  velocity  profile  is 
integrated  from  the  relation: 


^  =  1 

dy 


(18) 


The  velocity,  u,  is  the  local  velocity  direction  tangential  to  the  wall,  and  y  is  normal  to  the  wall. 
The  turbulent  viscosity  is  an  algebraic  function  of  the  local  shear  stress  and  the  length  scale. 

|I,  =  KypItl'^^Z)  (19) 


The  near  wall  damping  D,  is  a  function  of  a  friction  velocity  obtained  from  the  local  shear  stress, 

1  /  2 

Equation  18,  ii^  =  jx/pl  ,  and  the  damping  function  of  Equation  16. 

The  wall  function  is  evaluated  using  the  secant  approximation  to  Newton’s  method  to  solve 
for  V  The  iteration  continues  at  each  boundary  point  until  a  wall  shear  stress  found  that  gives  a 
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velocity  at  the  cell  center  that  matches  the  current  cell  center  velocity  in  the  finite  volume  solu¬ 
tion.  At  each  finite  volume  iteration  one  or  two  secant  iterations  are  normally  required  for  accept¬ 
able  convergence.  Each  secant  iteration  step  requires  numerical  integration  of  Equation  18,  which 
is  performed  on  a  one-dimensional  mesh  with  about  20  mesh  points.  Use  of  the  wall  function  iter¬ 
ation  reduces  the  near  wall  spacing  requirements  by  roughly  an  order  of  magnitude  in  an  LES 
solution.  In  addition,  a  coarser  near  wall  mesh  may  allow  larger  computational  time  steps.  How¬ 
ever,  because  the  wall  function  does  not  resolve  the  fine  scale  eddies  in  near  the  wall,  it  should  not 
be  used  where  a  high  fidelity  LES  boundary  layer  solution  is  required. 

2.7  Numerical  Diffusion  and  Limiting 

LES  predictions  are  extremely  sensitive  to  numerical  diffusion.  The  implicit  numerical  dif¬ 
fusion  in  lower  order  flux  schemes  can  be  of  the  same  order  as  the  sub  grid  scale  diffusion  from 
the  LES  equations.  The  weapons  bay  problem  is  particularly  challenging  for  LES  because  of  the 
strong  moving  shocks  responsible  for  the  acoustic  phenomena  in  the  bay.  Numerical  stability  for 
flows  with  strong  shocks  typically  requires  the  application  of  upwind  schemes  with  flux  limiters. 
Unfortunately,  the  upwind  schemes  typically  include  significant  numerical  diffusion.  For  this  pro¬ 
gram,  several  approaches  were  investigated  for  accurate  compressible  flow  simulations.  A  third 
order  scheme  was  implemented  in  parallel  Splitflow.  Two  hybrid  approaches  for  reducing  the 
implicit  diffusion  of  Roe’s  scheme  were  developed  and  tested.  Implementation  of  an  Essentially 
Non-Oscillatory  methocfwas  investigated  for  Splitflow 

Third  order  scheme:  A  third  order  accurate,  upwind  discretization  of  the  convective  flux  terms 
was  implemented  in  Splitflow.  This  scheme  includes  a  flux  limiter  that  makes  the  scheme  Total 
Variation  Diminishing  (TVD).  While  the  limiter  helps  to  prevent  spurious  oscillations,  particu¬ 
larly  in  flows  with  shock  waves,  it  adds  implicit  diffusion  to  the  solution.  Fluxes  are  constructed  at 
each  grid  face  in  the  flow  field  by  extrapolating  the  primitive  variables,  (p,  u,  v,  w,  T)  from  cell 
centers  to  cell  faces  to  obtain  right  and  left  extrapolated  values  at  each  face.  For  regions  of  the 
flowfield  where  the  mesh  is  regularly  spaced,  the  Splitflow  and  structured  mesh  results  are  identi¬ 
cal.  Notation  for  the  primitive  variables  and  faces  is  shown  in  Figure  2.15.  For  an  unlimited  flux, 
the  third  order  scheme  is: 
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Figure  2.15  Cell  notation  for  flux  descriptions.' 
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Limiting  is  included  in  this  scheme  with  a  compression  parameter,  c: 


(20) 


1 

'"2 

+ 


^i-\ 


+  ]^minmod[{q-_  ,  -  ^,  _2)>  c(^,-  -  _  i)] 

(9,  -^,  -  l)] 


(21) 


The  compression  factor  should  be  set  between  1  and  4.  Setting  the  factor  to  4  is  the  most  compres¬ 
sive  option  that  maintains  the  TVD  property,  resulting  in  the  smallest  truncation  error  in  the  solu¬ 
tion.  Setting  the  compression  factor  to  1  reduces  this  scheme  to  a  second  order  “minmod”  scheme. 
The  minmod  function  is  the  standard  defipition,  where 

a,  if{\a\  ^  1^1,  ab>0) 

minmod{a,  b)  =  b,if{\a.\<\b\,ab>Q)  (22) 

0,  if{ab  <  0) 

Construction  of  the  fluxes  on  the  i+1/2  face  is  similar.  Construction  of  the  fluxes  is  not  as  simple 
in  a  Splitflow  mesh  as  it  is  in  a  structured  grid  mesh  because  the  cell  face  for  one  cell  can  be  split 
into  multiple  faces  for  neighboring  cells,  as  shown  in  Figure  2.16.  For  a  one  dimensional  struc¬ 
tured  mesh,  evaluation  of  3rd  order  upwind  flux  at  a  face  requires  a  four  point  stencil  with  two 
cells  on  either  side  of  each  face.  Since  the  center  of  a  Splitflow  cell  is  no  longer  aligned  on  a  nor- 
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Figure  2.16  Primitive  variable  extrapolation  scheme  considers  tangential  components 
mal  to  the  face  center,  gradients  in  direction  parallel  to  cell  faces  must  be  accounted  for  in  flux 
extrapolation.  Cell  center  gradients  of  the  primitive  variables  are  evaluated  using  Gauss  theorem, 

V,  =  i  (23) 

where  A  is  the  outward  facing  normal  vector,  C2  is  the  cell  volume,  and  {qo-m)interp  •^he  value  of 
q  interpolated  between  the  cell  “0”  and  the  cell  m.  The  sum  of  the  dot  product  of  the  vector  r  and 
Vq  gives  a  second  order  accurate  value  of  the  primitive  variable  at  the  cell  face.  To  use  the 
upwind  formulation  of  Equation  21,  an  extrapolation  representative  of  the  difference  in  q  between 
cells  0  and  6.  and  an  interpolation  representative  of  the  difference  in  q  between  cells  0  and  3  must 
be  determined.  By  subtracting  the  contribution  of  the  interpolated  difference  between  cells  0  and 
3  from  the  cell  centered  gradient  term,  the  extrapolated  difference  is  derived.  Equation  21  is  then 
applied  to  these  differences  to  obtain  the  right  and  left  primitive  variables  at  the  face.  This  method 
gives  third  order  accuracy  for  regions  of  the  grid  with  equally  spaced  cells.  For  the  portions  of  the 
mesh  with  varying  spacing,  the  formal  accuracy  reverts  to  second  order.  The  truncation  error  in 
regions  of  the  grid  with  variations  in  the  cell  size  has  a  smaller  truncation  error  than  either  a  fully 
upwind  extrapolation  or  a  central  difference  extrapolation. 

Roe  Scheme  Limiting.  Fluxes  at  cell  faces  are  obtained  from  the  left  and  right  extrapolations  of 
the  primitive  variables. 
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F,{Q)  =  ^[F,(Qi^)  +  F,{Q,)-<i,\A\(Qi!-QL)]  <24) 

where  Ql  and  Qr  are  the  flow  variables  from  the  cells  to  the  left  and  the  right  of  the  face,  and 
A=  (3F/3j2).  The  Roe  averaged  state,  A  is  obtained  by  evaluating  A  using  Roe  averaged  val¬ 
ues  from  the  left  and  right  states.^  The  matrix  |a|  is  obtained  by  pre-  and  post-multiplying  the 
absolute  value  of  the  diagonal  matrix  of  eigenvalues  of  A  by  the  diagonalizing  matrices  obtained 
from  the  eigenvectors  of  A.  The  coefficient  0),  the  “Roe  Factor”,  controls  the  amount  of  Roe 
scheme  based  diffusion  to  be  included  in  the  simulation.  For  Roe  scheme,  the  factor  co  =  1,  while 
for  central  differencing,  (O  =  0.  In  the  “Large  Eddy  Simulation  for  Cavity  Flow  Analysis”  pro¬ 
gram,  the  Roe  flux  difference  splitting  technique  was  found  excessively  diffusive  for  a  low  speed, 
Mach  0.2,  channel  flow  simulation.  The  matrix  |a|  represents  a  matrix  of  smoothing  coefficients 
for  the  Euler  equations.  Because  some  of  the  coefficients  are  functions  of  the  speed  of  sound, 
these  coefficients  can  dominate  the  diffusion  for  low  speed  flow.  The  goal  for  the  current  program 
was  to  assess  the  effectiveness  of  Roe  scheme  for  a  higher  Mach  number  flow,  and  modify  Roe 
scheme  as  required  for  transonic  and  supersonic  applications. 

Higher  speed  flow  simulations  should  be  less  sensitive  to  Roe  scheme  because  the  damping 
coefficients  in  the  Roe  matrix  proportional  to  the  speed  of  sound  are  similar  in  magnitude  to  the 
flow  velocity  in  this  case.  A  Mach  1.2  channel  flow  simulation  was  performed  for  this  program. 
From  multiple  tests,  it  was  found  that  while  some  unsteady  turbulent  structures  could  be  obtained 
with  CO  =  I,  higher,  more  realistic  levels  of  resolved  turbulence  were  obtained  for  to  =  0.2.  In  r.ddi- 
tion,  use  of  the  flux  limiter  proved  necessary  to  obtaining  a  stable  flow.  The  compression  factor 
was  .set  to  4.0  for  the  channel  flow  simulation,  the  most  compressive  .setting.  In  contrast,  the  low 
speed  solution  required  0)  =  0.1  to  obtain  reasonable  turbulence  levels,  and  the  flux  limiter  was 
turned  off.  Table  2. 1  shows  the  differences  in  the  settings  between  the  low  speed  solution  obtained 
with  the  hybrid  version  of  Splitflow  and  the  higher  speed  solution  obtained  with  parallel  Splitflow. 

Table  2.1  Flux  scheme  coefficient  settings  for  low  speed  and  transonic  simulations 


Simulation 

Mach 

Roe 

Number 

Factor 

Limiting  Compression  Factor 

0.2 

0.1 

No  limiting 

1.2 

0.2 

4.0 
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As  expected,  the  higher  speed  solutions  could  tolerate,  and  in  fact  required,  larger  Roe  coefficients 
and  flux  limiting  than  the  low  speed  solutions. 

An  alternative  approach  for  Roe  Scheme  limiting  was  developed.  In  this  scheme,  instead  of 
limiting  the  entire  Roe  Scheme  diffusion  matrix,  only  one  speed  of  sound  related  term  was  lim¬ 
ited.  The  product  of  the  Roe  diffusion  matrix  is  and  the  primitive  variable  differences, 

[ap  Am  Av  Aw  Ap] 


1 

0 

u 

Am  -  q^A(/^ 

|Af,  3.4I  = 

^  a  ^ 

V 

w 

-2  -2  -2 
M  +  V  +  W 

+  P 

Av  -  q^Af/^ 

Aw  -  q.At/^ 

mAm  +  vAv  +  wAw  -  UcAU^ 

L 

2 

(25) 


1 


Uc±a 


[Ap  ±  pa(ri^AM  +  rj^.Av  +  ri,Aw)] 
2a^ 


u  ±  pq 

V  ± 

w  ±  pq.M 

h,  ±  U 


where  U ^  =  ur\  ^  +  vq^,  +  wq.  is  the  contravarient  component  of  velocity.  For  the  standard  Roe 
.scheme,  P  =  1 .  The  term  for  which  P  is  a  coefficient  is  proportional  to  the  speed  of  sound.  By  lim¬ 
iting  this  coefficient.  Roe’s  scheme  can  be  less  diffusive  for  LES  while  reducing  spurious  flow 
o.scillaiions  at  shock-waves. 

Roe's  scheme  is  based  on  the  Euler  equations  in  one  dimension.  When  the  flow  direction  is 
close  to.  but  not  exactly  parallel  with  a  cell  face  and  a  shear  strain  is  present  in  the  direction  nor¬ 
mal  to  the  face,  the  [AFj  5]  term  of  the  Roe  diffusion  vector  can  have  a  significant  effect  on  the 
shear  stress.  For  LES  simulations,  this  term  adds  excessive  damping  to  shear  and  vortical  regions 
of  the  flow  domain.  Conversely,  when  the  flow  direction  is  nearly  normal  to  the  cell  face,  differ¬ 
ences  in  the  velocity  across  the  face  is  primarily  an  inviscid  phenomenon  and  the  Roe  scheme 
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based  on  the  Euler  equations  is  appropriate.  For  the  modified  scheme,  p  is  set  to  a  function  that 
reduces  the  diffusion  magnitude  when  the  flow  direction  is  not  aligned  with  the  face  normal. 


P 


=  maximum 


\U. 


o. 


/  2  ,  2  ^  2  , 
VM  +  V  +  W  J 


(26) 


The  constant  a  was  found  to  be  necessary  to  maintain  some  of  the  effects  of  this  smoothing  coef¬ 
ficient  even  when  the  flow  nearly  parallel  to  the  cell  face.  Typically,  we  set  a  =  0.1.  Results 
obtained  with  this  model  will  be  discussed  in  Section  3. 

ENO  scheme  evaluation:  Recently,  Essentially  Non-Oscillatory  (ENO)  and  Weighted  ENO 
(WENO)  schemes^  have  been  used  for  finite  difference  and  finite  volume  LES  applications.  ENO 
schemes  provide  a  procedure  for  obtaining  high  order  accuracy  for  piecewise  smooth  problems 
with  discontinuities.  An  ENO  scheme  could  be  used  in  place  of  the  extrapolation  and  limiting  pro¬ 
cedures  specified  in  Equations  20,  21,  and  22.  A  Roe  flux  difference  split  method  can  be  used  to 
determine  the  flux.  The  ENO  procedure  can  provide  higher  order  expressions  for  unevenly  spaced 
structured  grids.  As  part  of  this  program,  the  implementation  of  a  third  order  ENO  scheme  in 
Splitflow  was  investigated. 

A  third  order,  upwind  ENO  scheme  requires  information  from  three  cells  on  each  side  of  a  flux 
face.  By  comparison,  the  Splitflow  third  order  scheme  requires  two  cells  on  each  side  of  a  face. 
For  an  unstructured  grid  with  split  cells  such  as  Splitflow,  expansion  of  the  flux  stencil  adds 
greatly  to  the  computational  complexity  of  the  scheme.  A  method  to  implement  an  ENO  scheme 
in  Splitflow  was  developed.  This  scheme  uses  the  cell  center  gradients  at  two  cells  on  either  side 
of  the  cell  face  to  construct  the  ENO  differences.  In  the  logical  structure  of  Splitflow  only  cells 
touching  a  cell  face  can  be  readily  identified.  For  the  current  third  order  scheme,  flux  evaluation  is 
a  two  step  process,  where  cell  center  gradients  are  determined  using  one  pass  through  all  cell 
faces,  and  the  extrapolations  and  fluxes  are  obtained  by  a  second  pass  through  all  cell  faces  in  the 
flow  domain.  While  a  method  for  identification  of  the  additional  cell  center  information  required 
for  ENO  was  developed,  the  method  was  unwieldy  for  parallel  implementation.  The  difficulties 
arose  from  the  need  to  pass  excessive  amounts  of  information  large  numbers  of  times  across 
boundaries  between  the  flow  domains  on  each  Splitflow  parallel  node.  A  practical  method  for 
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implementation  of  and  ENO  scheme  in  Splitflow  would  require  a  total  rewrite  of  the  Splitflow 
data  structure.  This  is  a  task  that  would  require  many  man  months  to  complete,  so  this  approach  to 
improve  numerical  accuracy  was  not  pursued  further.  However,  it  does  appear  that  implementa¬ 
tion  of  an  ENO  scheme  is  feasible,  and  results  from  the  literature  indicate  that  implementation  of 
an  ENO  scheme  could  benefit  accuracy  and  stability  of  Splitflow  for  LES  applications. 

2.8  Periodic  Boundary  Conditions 

For  channel  and  cavity  flow  studies,  periodic  boundary  conditions  allow  small  computa¬ 
tional  domains  to  be  used  for  the  efficient  development  of  turbulent  boundary  layers.  Periodic 
boundary  conditions  allow  the  turbulent  structures  to  be  passed  across  the  boundaries  of  the  flow 
domain.  Two  types  of  periodic  boundaries  were  required  for  this  research,  a  lateral  periodic 
boundary  condition  and  an  inflow  -  outflow  periodic  condition.  The  inflow  -  outflow  periodic 
boundary  condition  used  in  this  research  underwent  significant  modifications  compared  to  the 
boundary  condition  used  in  the  “LES  for  Cavity  Flow  Analysis”  program.  In  the  earlier  research, 
the  boundary  condition  was  applied  to  a  Mach  0.2  channel  flow.  This  formulation  did  not  work 
well  for  the  Mach  1.2  channel  flow  in  the  current  research.  The  new  formulation  worked  well. 

Periodic  Boundary  Condition  Flux  Calculation 

For  cut  grids,  periodic  boundaries  are  defined  by  the  triangular  “geometry  facets”.  The 
geometry  facets  on  each  side  of  the  periodic  boundary  should  match  exactly  in  size,  shape  and  ori¬ 
entation.  The  geometry  facets  are  cut  oi.  each  side  into  “boundary  facets”  by  the  Cartesian  cell 
faces.  Flow  variables  are  accumulated  using  an  area  weighted  procedure  from  the  boundary  facets 
to  their  corresponding  geometry  facets  from  each  side  of  the  flow  domain.  Fluxes  are  then  calcu¬ 
lated  on  the  geometry  facets.  The  fluxes  from  the  geometry  facets  are  then  distributed  to  the 
appropriate  cells  using  the  boundary  facets.  Figure  2.17  shows  the  relationship  between  geometry 
facets,  cell  faces  and  boundary  facets.  The  accuracy  of  the  periodic  boundary  condition  depends 
on  the  resolution  of  the  geometry  facets. 

For  body  fitted  Omnigrids,  geometry  facets  are  generated  by  triangulating  the  quadrilateral 
faces  on  the  periodic  boundary.  As  a  result,  each  cell  face  is  associated  with  two  geometry  facets 
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Figure  2.17  Periodic  Boundary  Facet  Definitions 

that  fit  perfectly  within  its  face.  As  a  result,  the  periodic  boundary  condition  is  more  accurate  for 
Omnigrid  calculations. 


Compressible  Inflow  •  Outflow  Periodic  Boundary  Conditions 

For  fully  developed  channel  flow,  jump  conditions  must  be  included  between  the  periodic 
inflow  and  outflow  boundaries.  We  hold  mass  flux  and  total  energy  flux  constant  at  the  inflow,  and 
the  average  static  pressure  constant  at  the  outflow.  The  flow  variables  from  the  outflow  must  be 
adjusted  to  be  used  at  the  inflow  and  the  inflow  must  be  adjusted  for  use  at  the  outflow.  A  sche¬ 
matic  that  shows  the  notation  for  the  periodic  boundary  condition  derivation  is  shown  in 
Figure  2.18.  Surface  integrals  of  inflow  and  outflow  pressure  and  outflow  mass  flux  and  total 
energy  flux  are  evaluated.  The  average  pressure  is: 

^0  =  ^0  h  -  dz  (27) 
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Figure  2.18  Diagram  of  notation  for  periodic  inflow/outflow  boundary  condition. 
The  ma.ss  flux  integrals  are  defined  as: 


h  =  Ji  =  \\p2i^2^y 

The  total  energy  integrals  are  defined  as 

^0  =  ^2  =  Jjp2«2^r2^>'^-  (29) 

\  1  2 

The  total  energy,  Ej,  is  the  sum  of  the  internal  energy,  e,  and  -(«“  +  v'  +  vv-  ) .  For  the  LES  for 
Weapons  Bay  Analysis  program,  total  enthalpy  was  used  instead  of  total  energy.  For  the  inflow 
conditions  at  each  geometry  facet,  assuming  the  perfect  gas  equation  of  state,  we  have  five 
unknowns.  Pi/,,  density,  r/y^,  vyy,,  and  wyy,,  velocity  components,  and  e^,,  so  five  equations  in 
terms  of  the  outflow  variables  and  integrals  are  required.  We  a.ssume  Vjy,  =  V2(M,y,/ Uj)  and 
U  =  us(//i/,/»t)  .  Then  from  the  pressure,  mass  and  total  enthalpy  equations,  we  obtain  three 
additional  equations. 


P\b  ~  ^2  ^0  “  ^2 


^*\h[P\h^ \h'^  P \h'^  ^P\b^ \h\  -  “2(  P2^2  +  ^2  7P2^2 


1  . 


(30) 
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where  and  V2  =  m?  +  v?  +  w]  .  Substituting  the  first  two  equations 

from  Equation  30  into  the  third  equation,  a  quadratic  equation  for  m;/,  is  obtained 


P2V27, 


'2^14. J  j 


0  2 

^\b 


+  h^^\b~  P2"2^r2^ 


^0 

=  0 


(31) 


After  solving  for  u we  can  then  solve  for  the  other  variables. 

A  similar  approach  is  used  for  the  outflow  periodic  boundary.  We  assume 
V2^  =  V|(t/7^/w,)  and  w,^  =  w,(M2i,/Mi),  and  obtain  for  the  pressure,  mass  and  total 
enthalpy  equations: 


^2b  -  ^1  ^2  “  ■^l 

Plb'^lb  -  Pl'^1  / 

•^0 

^*2h{p2b^2b'^  ^2b'^  2^2b^2hj  ~  ^  I  (^P  1  ^  I  ^1  2^  1 


(32) 


As  for  the  inflow  case,  we  now  obtain  a  quadratic  equation  for  «2/>- 


P  1  V 1  J ')  7  V  ^  ^ 

~^^Y‘*~2b'''  y — r(^i +^2"^i)"2/j”Pi"i^nl^  =  ® 
2»i  -/o  Y -  *  ^0 


(33) 


The  pressure  relations  for  and  in  Equations  30  and  32  are  different  than  the  relation  used 
in  the  LES  for  Cavity  Flow  Analysis  program.  The  old  method  tended  to  over-specify  the  pres¬ 
sure.  The  inflow  -  outflow  periodic  boundary  condition  succeeded  in  passing  turbulent  flow  struc¬ 
tures  with  limited  distortion  from  the  outflow  into  the  inflow  of  the  channel  flows. 
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Section  3:  Demonstration  Cases 

Three  Splitflow  demonstration  cases  were  run  as  part  of  this  program.  A  turbulent  channel 
flow  is  used  to  provide  the  inflow  turbulent  boundary  layer  to  a  simple  cavity  geometry.  An  F-22 
Store  Release  simulation  was  also  performed.  These  cases  represent  first  attempts  at  calculating 
these  complex,  compressible  flows  with  Splitflow.  The  accuracy  and  computational  robustness  of 
these  demonstration!:;  was  not  as  good  as  we  had  hoped  they  would  be.  In  particular,  the  turbulent 
boundary  layer  results  from  the  channel  simulation  are  degraded  by  excessive  numerical  diffu¬ 
sion.  The  parallel  scalability  of  the  automatic  gridding  logic  is  worse  than  that  of  the  solver,  and 
its  frequent  use  to  update  the  geometry  caused  the  overall  scalability  to  degrade  on  large  numbers 
of  processors.  As  a  result,  simulations  took  excessive  wall  clock  time.  In  addition,  as  the  store 
moved  through  the  grid,  the  grid  cutting  process  occasionally  resulted  in  problem  cells  which 
caused  numerical  instability  problems.  Due  to  these  computational  problems,  the  duration  of  the 
F-22  store  release  simulation  was  limited  and  the  final  position  of  the  store  in  the  simulation  is  not 
tremendously  far  from  the  initial  position.  While  disappointing,  these  results  provide  valuable 
information  on  both  the  strengths  and  weaknesses  of  the  current  procedures  and  where  additional 
development  work  should  be  focused. 

3.1  Channel  Flow 

The  channel  flow  provides  a  good  test  case  for  near  wall  compressible  LES  capability.  In 
addition,  the  channel  flow  case  was  used  to  provide  an  inflow  boundary  layer  for  the  simple  cavity 
flow  solution.  The  Mach  number  at  the  channel  centerline  1.2.  The  physical  domain  is  one  inch 
wide,  two  inches  long  and  0.56  inches  high.  The  channel  height  is  twice  the  thickness  of  the 
incoming  boundary  layer  for  the  cavity  flow  problem.  The  Reynolds  number  for  the  channel  using 
one  half  the  channel  height  as  the  characteri.stic  length  is  Re^  =  19000. 

The  grid  was  generated  using  the  cut  grid  option  in  Splitflow.  The  near  wall  grid  spacing  is 
wall  normal,  Az"*"  =  1.7,  lateral,  Ay'^  =  29  and  streamwise,  Ax"^  =  58.  A  side  view  of  the  grid  is 
shown  in  Figure  3.1,  while  a  cross  section  view  is  shown  in  Figure  3.2.  Figure  3.3  is  a  zoomed-in 
view  of  the  grid  in  the  cross  flow  direction  near  the  wall  showing  the  additional  levels  of  cell 
refinement  employed  near  the  wall.  Similar  refinement  is  employed  in  the  streamwise  direction 
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near  the  wall.  The  total  number  of  cells  for  the  simulation  is  138,240.  If  the  grid  had  been  gener¬ 
ated  with  a  structured  mesh  with  the  same  near  wall  spacing  would  require  327,680  cells. 

Results  using  the  two  different  types  of  Roe  scheme  limiting,  described  in  Section  2.7,  will 
be  presented.  For  the  first  case,  the  Roe  scheme  factor  (O  =  0.2,  while  for  the  second  case,  the 
modified  Roe  scheme  based  on  Equation  26  is  used  with  a  =  0.1.  The  channel  flow  velocities. 


Figure  3.4  Channel  flow  velocity  profiles  for  two  versions  of  Roe  scl^eme. 


shown  in  Figure  3.4  appear  reasonable.  However,  the  profiles  do  not  agree  well  with  the  law  of  the 
wall.  Velocity  profiles  from  the  two  schemes  are  eompared  in  Figure  3.5.  The  modified  Roe 
scheme  gives  slightly  more  accurate  results  than  simply  limiting  the  Roe  scheme  diffusion  with 
the  Roe  factor.  Profiles  were  evaluated  using  both  incompressible  and  compressible  wall  units.  A 
comparison  of  incompressible  and  compressible  wall  units  for  the  modified  Roe  scheme  case  is 
shown  in  Figure  3.6.  Compressible  wall  units  account  for  the  effect  of  density  variations  in  the 
boundary  layer  using  the  van  Driest  transformation. 

The  resolved  turbulent  energy  was  evaluated  for  the  channel  flow  calculation.  The  turbulent 
intensities  are  shown  in  Figure  3.7.  The  peak  level  for  the  streamwise  turbulent  intensity  is  signif¬ 
icantly  lower  than  the  peak  turbulent  intensity  levels  found  experimentally.  Based  on  the  data  of 
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Wei  and  Willmarth,^  the  ratio  of  the  streamwise  peak  intesity  to  the  friction  velocity  is  approxi¬ 
mately  2.5,  and  the  ratio  of  intensity  normal  to  the  wall  to  the  friction  velocity  is  approximately 
1.0.  Our  resolved  intensities  are  u'lu^  =  1.4  and  w'lu^  =  0.15. 


Figure  3.7  Turbulent  intensities  for  channel  flow 
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3.2  Cavity  Flow 

The  cavity  flow  demonstration  is  useful  for  assessing  the  effectiveness  of  Splitflow  with  the 
LES  model  for  acoustic  load  simulation.  A  flow  case  from  the  Weapons  Internal  Cavity  and  Sepa¬ 
ration  (WICS)  program  was  chosen  for  comparison.^  The  test  model  has  a  sharp  leading  edge  and 
a  15  inch  long,  16.5  inch  wide  flat  plate  in  front  of  the  cavity.  The  cavity  L  x  W  x  D  =  18”  x  4”  x 
4”.  The  flow  condjt-.ons  for  the  case  selected  are  summarized  in  Table  3.1. 

Table  3.1  Flow  conditions  for  cavity  simulation 
Mach  1.2 

U„  1208  ft/sec 

T„  422  R 

Poo  188  psf 

qoo  189  psf 

1.0x10^ 

In  order  to  obtain  a  good  quality  incoming  turbulent  boundary  layer  that  includes  unsteady 
structures,  a  small  channel  domain  was  defined  and  run  with  periodic  boundary  conditions.  The 
methods  for  this  simulation  are  described  in  Section  2.4  and  Section  2.8  and  the  results  are  dis¬ 
cussed  in  Section  3.1.  The  one  inch  wide  channel  domain  was  applied  across  the  computational 
domain  by  duplicating  the  solution  along  the  inflow.  Because  the  channel  flow  solution  uses  peri¬ 
odic  boundary  condition  in  the  lateral  direction,  the  flow  variables  are  continuous  along  the  inflow 
plane.  The  distanc"  between  the  center  of  the  channel  and  the  floor  of  the  channel  corresponds  to 
the  thickness  of  the  incoming  boundary  layer. 

The  cavity  mesh  was  generated  using  the  omnigrid  mesh  generator.  The  cavity  grid  con¬ 
tained  614,000  cells.  A  centerline  view  of  the  mesh  is  shown  in  Figure  3.8.  A  cross  section  cut 
through  the  grid  two  inches  downstream  of  the  cavity  leading  edge  is  shown  in  Figure  3.9. 

As  expected,  a  highly  unsteady  flowfield  was  obtained  in  the  cavity.  In  addition,  significant 
asymmetry  was  observed.  Contours  of  the  streamwise  velocity  component  at  one  instant  of  time 
in  the  solution  along  the  cavity  centerline  are  shown  in  Figure  3.10.  Note  the  significant  region  of 
strong  reverse  flow  present  in  the  solution.  Streamwise  velocity  contours  at  cross  section  cuts  2, 
10  and  16  inches  down.stream  of  the  cavity  leading  edge  are  shown  in  Figures  3.1 1,  3.12  and 
3.13. 
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Figure  3.8  Centerline  view  of  cavity  grid  generated  using  omnigrid. 


The  flow  solution  was  computationally  intensive.  Stable,  accurate  solution  of  the  incoming 
channel  flow  required  use  of  a  time  step  cl  5x10"^  second.  This  is  not  unreasonably  small  for  the 
channel  flow  solution,  since  it  represents  a  CFL  =  (U  +  a)At/(Ax)  =  0.4  in  the  streamwise 
direction.  The  CFL  number  in  the  direction  normal  to  the  wall,  due  to  the  speed  of  sound  compo¬ 
nent  a  in  the  CFL  definition  and  the  fine  grid  in  the  wall  normal  direction,  is  several  orders  of 
magnitude  larger.  As  a  result,  many  thousands  of  iterations  were  required  to  obtain  a  Power  Spec¬ 
tral  Density  (PSD)  spectra.  Power  Spectral  Density  results  were  generated  for  four  of  the  loca¬ 
tions  where  high  response  pressure  data  were  available.  The  spectra  were  based  on  16384  samples 
averaged  over  two  samples.  This  results  in  a  bandwidth  of  121  Hz.  As  expected,  the  highest 
acoustic  levels  are  found  near  the  back  wall  of  the  cavity.  Cavity  geometry  and  location  of  the 
pressure  transducers  are  given  in  Figure  3.14.  Acoustic  spectra  are  shown  in  Figures  3.15,  3.16, 
3.17  and  3.18. 
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Figure  3.9  Cross  section  cut  through  the  cavity  grid. 
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Figure  3.10  Velocity  contours  on  cavity  centerline. 
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Figure  3.12  Velocity  contours  10  inches  from  leading  edge  of  cavity 


Figure  3.13  Velocity  contours  16  inches  from  leading  edge  of  cavity,  2  inches  upstream  of  back  wall. 
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Figure  3.16  Acoustic  Spectra  for  KIO,  5.575”  from  front  wall  of  cavity 


Frequency  (Hz) 

Figure  3.17  Acoustic  Spectra  for  K13, 10.075”  from  front  wall  of  cavity 
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3.3  F-22  Store  Release  Simulation 


Release  Conditions 


A  release  condition  was  chosen  for  this  simulation  from  the  F-22  program  wind  tunnel  data 
base.  The  chosen  case  is  interesting  and  challenging  since  the  store  shows  significant  response  to 
the  aircraft  flow  field.  Configuration,  flow  and  release  conditions  for  the  case  are  described  in 
Table  3.2.  The  flow  field  and  trajectory  were  simulated  1/15  scale  (model  scale).  Note  that  the 
load  factor  (Nz)  is  less  than  1.0  meaning  that  a  push-over  maneuver  was  modeled. 


Table  3.2'najectory  simulation  release  conditions 


Variable 

Value 

Test 

4T-2 

Run# 

1756 

Mach 

1.6 

Altitude 

20.000  Ft 

Nz 

+0.5  g 

a 

1  deg 

b 

Odeg 

'^EOS 

28  Ft/sec 

Qeos 

-0.5  rad/sec 

Station  7 

Empty 

Station  8 

Empty 

Station  9 

AIM-120C 

Grid  Generation 

Two  grid  domains  were  used  for  the  store  release  simulation.  A  body  conforming  omnigrid 
was  generated  at  the  forebody  of  the  aircraft  in  order  to  obtain  accurate  boundary  layer  profiles  at 
the  weapons  bay  leading  edge.  A  cut  grid  surrounded  this  grid  and  the  remainder  of  the  aircraft. 
More  details  on  the  grid  interface  are  contained  in  Section  2.4,  Interface  Grid  Methods.  A  view  of 
the  F-22  geometry  that  was  resolved  and  simulated  is  shown  in  Figure  3.19. 
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Figure  3.19  Full  view  of  F-22  geometry  simulated  and  Aim  120  store  in  initial  position. 

Initial  Conditions 

The  EOS  store  CG  location  was  9  vertical  inches  (+Z  in  body  axes)  from  the  full  scale  car¬ 
riage  position.  The  store  orientation  was  -1.67  degrees  nose  down  in  the  body  axis  system.  This 
resulted  in  an  initial  position  of  (533.787,  30.415,  59.210)  inches  full  scale  and  (2.965,  0.1689, 
0.3289)  feet  model  scale.  The  store  CG  at  the  carriage  position  is  (533.787,  30.415,  68.210)  full 
scale  inches  and  (2.965,  0.1689,  0.3789)  feet  model  scale.  The  EOS  linear  and  rotational  veloci¬ 
ties  (Table  3.2)  were  applied  as  an  initial  condition  for  the  store. 

Moving  Body  Solutions 

Figure  3.20  shows  the  grid  in  the  region  near  the  interface  between  the  boundary  fitted  grid 
and  the  cut  grid  and  in  the  region  of  the  moving  store  for  different  time  levels.  The  addition  of  grid 
cells  is  clearly  visible  as  the  store  moves  down  and  leaves  behind  a  trail  of  small  grid  cells  which 
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were  required  to  resolve  previous  geometry  positions.  Occasional  cell  deletion  was  programmed 
to  eliminate  cells  as  their  number  grew  excessively.  Figure  3.21  shows  another  view  of  the  simula¬ 
tion  and  illustrates  how  the  pressure  coefficients  on  the  store  are  changing  as  it  encounters  higher 
speed  flow  from  the  freestream  as  it  moves  downward  in  4.76  milliseconds  (model  scale).  The 
interface  between  the  grids  can  be  seen  below  the  upstream  edge  of  the  weapons  bay. 


Figure  3.10  Side  views  of  moving  store.  Both  grids  visible. 

Trajectory  Results 

Although  a  full  near-field  trajectory  was  not  obtained  due  to  computational  issues  discussed 
above,  initial  results  look  quite  good.  Figure  3.22  shows  sequences  of  the  first  0.00476  seconds 
model  scale  of  the  trajectory  (0.0714  seconds  full  scale),  with  flow  field  features.  An  area  of  low 
Mach  number  flow  is  observed  above  the  missile  revealing  an  influence  of  body  velocity  on  the 
flow  field.  A  plot  of  position  versus  time  is  shown  in  Figure  3.23.  Position  is  in  full  scale  feet  and 
time  is  in  model  scale  seconds  for  comparison  to  the  above  flow  field  figures. 
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Figure  3.21  Two  store  positions  indicate  changes  in  aerodynamic  forces. 
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Mach  Number  Pressure  Coefficient 


Figure  3.22  Store  trajectory  with  mach  and  pressure  contours  in  flow  field. 
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Section  4:  Conclusions 

Significant  progress  was  made  toward  the  development  of  a  highly  accurate  method  that 
integrates  acoustic  load  prediction  and  store  separation  simulation  using  CFD.  The  program  goals 
were  extremely  ambitious.  While  more  work  is  needed  improve  the  computational  efficiency  of 
this  highly  automated  method,  this  program  pioneered  the  development  and  application  of  this 
methodology.  The  main  accomplishments  of  this  program  were: 

1 .  Developed  moving  body  capability  in  the  highly  automated,  Splitflow  parallel 
unstructured  grid  flow  solver 

2.  Integrated  six  degree  of  freedom  dynamics  capability  into  the  Splitflow  solver 

3.  Developed  a  method  to  simultaneously  solve  multiple  Splitflow  solutions  separated  by 
faceted  surface  boundaries  that  pass  fluxes  conservatively  between  the  solutions 

4.  Demonstrated  the  capability  to  simulate  the  release  of  a  real  store  from  an  F-22 

5.  Predicted  the  unsteady  acoustic  loads  and  the  PSD  using  Splitflow 

6.  Developed  a  methodology  for  simulation  of  compressible  channel  flows  with  periodic 
inflow-outflow  boundary  conditions 

Much  was  learned  during  the  methodology  development  and  application  efforts  for  this  pro¬ 
gram.  These  lessons  point  to  areas  were  additional  work  is  required.  The  most  important  lessons 
from  this  study  are: 

1 .  Development  of  flux  discretization  schemes  that  are  high  order  accurate,  low 
diffusion  for  compressible  flow  analysis  using  unstructured  parallel  flow  solvers 
remains  a  significant  challenge.  High  accuracy  is  essential  for  accurate  Large  Eddy 
Simulation. 

2.  Improvements  to  the  parallel  scalability  of  Splitflow  are  needed  for  practical  application 
of  viscous  unsteady  flow  problems  around  complex  geometries  and  moving  bodies. 

3.  Numerical  approaches  for  initializing  the  flowfield  for  cells  created  at  the  boundaries  of 
a  moving  surface  undergoing  large  displacements  require  additional  development  to 
assure  efficient  and  accurate  results. 

The  capabilities  developed  in  this  research  effort  are  important  to  the  overall  goal  of  reduc¬ 
ing  costs  and  design  cycle  times  for  aerospace  programs.  Additional  development  effort  will  be 
required  before  these  methods  will  be  ready  for  routine  application.  However,  when  considered  in 
light  of  continued  improvements  in  computer  performance,  reductions  in  aircraft  program  devel¬ 
opment  budgets,  and  high  cost  of  wind  tunnel  testing,  the  need  for  weapons  bay  simulation  meth¬ 
ods  will  persist,  and  the  feasibility  of  their  use  will  continue  to  improve. 
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